On Approximation of the Eigenvalues of 
Perturbed Periodic Schrodinger Operators 

Lyonell Boulton and Michael Levitin 

Maxwell Institute for Mathematical Sciences 
and Department of Mathematics 
Heriot-Watt University 
Riccarton, Edinburgh EH14 4AS, U. K. 
email {L. Boulton, M.LevitinjOma. hw.ac.uk 
www.ma.hw.ac.uk/~lyonell/, www.ma.hw.ac.uk/~levitin/ 

February 2007 



Abstract 

This paper addresses the problem of computing the eigenvalues lying 
in the gaps of the essential spectrum of a periodic Schrodinger operator 
perturbed by a fast decreasing potential. We use a recently developed 
technique, the so called quadratic projection method, in order to achieve 
convergence free from spectral pollution. We describe the theoretical 
foundations of the method in detail, and illustrate its effectiveness by 
several examples. 



1 Introduction 

It is well known that the problem of approximating the eigenvalues lying in 
gaps of the essential spectrum of a self-adjoint operator by a sequence of finite- 
dimensional problems (e.g. for numerical analysis) is far from trivial. The pres- 
ence of essential spectrum both above and below an eigenvalue means that 
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there is no obvious variational principle (cf. e.g. [DoEsSe]), so an approxima- 
tion/computation by a standard projection method is not always possible. The 
main difficulty is due to the existence of sequences of eigenvalues of the (finite- 
dimensional) approximate operators, accumulating at points in the gaps which 
do not belong to the spectrum. These points are called spurious eigenvalues, 
and the phenomenon itself is often referred to as spectral pollution. 

It has been shown, for general unbounded self-adjoint operators, that spectral 
pollution in a projection method may occur at any real point of the resolvent set 
located between two parts of the essential spectrum (see [LeSh, Theorem 2.1]). 
This is a consequence of the fact that the resolvent is not compact. A substantial 
amount of research has been devoted to finding ways of choosing the projectors, 
in order to achieve a "safe" method for particular problems, see e.g. [RaSa^Va] 
and [BoBr]. Techniques vary considerably according to the problem and are by 
no means universal. 

In this paper we address the question of spectral pollution and its avoidance 
for a perturbed periodic Schrodinger operator 

H:=-A + V, (1.1) 

acting in the Hilbert space L^(M^), where V = Vp + Vd, with Vp being purely 
periodic with respect to some lattice of and Vd being fast decaying at infinity. 
The essential spectrum of H is determined by Vp. It consists of bands of abso- 
lutely continuous spectrum, separated by gaps in the resolvent set. If Vd = 0, 
the spectrum is purely essential. If Vd ^ 0, discrete eigenvalues may appear in 
the gaps, see [DeHe]. 

A usual method for finding the essential spectrum of H analytically, the so 
called Floquet-Bloch technique, has been well studied (see e.g. [ReSi], [Ku] and 
the references therein). It gives a decomposition of the periodic part of the 
operator into a direct integral of operators on a basic periodic cell. This reduces 
the problem of finding the endpoints of the bands in the essential spectrum, 
to the problem of finding the eigenvalues of differential operators in a compact 
domain with regular boundary conditions. 

Much less in known about the discrete spectrum of H, which has to be 
either estimated numerically or studied by means of asymptotic techniques (for 
the latter see e.g. [DeHe], [Bi] and [Su]). As we shall see below, the natural 
approach of truncating to a large compact domain and applying the projection 
method to the corresponding Dirichlet problem, is prone to spectral pollution. 
This makes the numerical localisation of these eigenvalues particularly difficult. 
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The purpose of this paper is to describe an alternative procedure for finding 
eigenvalues, the so called quadratic projection method, recently studied in an 
abstract setting in [Sh], [LeSh], [Bol] and [Bo2]. The distinctive feature of our 
method is that the underlying discretised eigenvalue problem is quadratic in the 
spectral parameter (rather than linear), and has non-real eigenvalues. Its main 
advantage over a standard projection method lies in its robustness: it never 
pollutes and it always provides a posteriori two-sided estimates of the error of 
computed eigenvalues. 

The paper is organised as follows. In Section 2 we discuss the phenomenon 
of spectral pollution in a standard projection method and discuss the quadratic 
projection in an abstract context. Our Corollary 2.6 is an improvement upon 
previously known non-pollution results for the general quadratic method. In Sec- 
tion 3 we provide details on how to implement the quadratic projection method 
for the numerical localisation of the eigenvalues of operator H. We also discuss 
some concrete numerical examples, but deliberately avoid including the full ac- 
count of the numerical procedures we have used, in order not to overload the 
text with unnecessary technical details. These will appear elsewhere. 

2 The quadratic projection method 

2.1 Spectral pollution in an ordinary projection method 

Before proceeding to describe our method, we want to give a rigorous motivation 
why it is needed at all, and why spectral pollution is intrinsic in the standard 
projection method. 

Let A be a self-adjoint operator in a Hilbert space H with a dense domain, 
Dom(A). The spectrum of A, Spec(A), can be decomposed into the discrete 
spectrum, SpeC(jig(.(74), consisting of isolated eigenvalues of finite multiplicity, 
and the essential spectrum, SpeCgss(A) := Spec(74) \ SpeCdisp(>l) = {A : ^4 — 
XI is not Fredholm}. 

Take a finite-dimensional subspace C C Dom{A), and let : H — > C be 
the orthogonal projection onto C Let Ac :— H^A \ C 

The projection method, also known as the Galerkin method, consists in trun- 
cating the (infinite-dimensional) spectral problem Au — \u to 

Acu = Alt for some u e C \ {0}. (2.1) 

If the operator A is bounded from below and has a compact resolvent, this 
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provides an effective way of estimating numerically the eigenvalues of A. The 
k-th eigenvalue of (2.1) will always be above the k-th eigenvalue of A, count- 
ing multiplicity, [ReSi, Section XIII. 1]. Furthermore, if £ approximates Dom(A) 
reasonably well, then the first few eigenvalues of (2.1) will be close to the corre- 
sponding ones of A. 

A precise statement can be easily obtained from the minimax principle: 

Lemma 2.1. Let Cn be a sequence of finite-dimensional subspaces ofT>om.[A). 
Assume that A is bounded below and has a compact resolvent. Let Ai < . . . < 
be the first m eigenvalues of A. Let 

£ — Span {it e Dom(^) : Au — XkU, 1 < k < m}, 

be the spectral subspace associated with {Ai, . . . , A^}. If 

lim \\AP{u-Uc^u)\\^0, (2.2) 

n— »oo 

holds for p = 0, 1 and all u E 8, then the k-th eigenvalue of (2.1) approaches 
the k-th eigenvalue of A as n —> oo for 1 < k < m. 

We omit the proof. 

In some particular cases it is also possible to estimate the convergence rate 
of the eigenvalues [StFi]. 

Similar results can be established if the resolvent of A is non-compact, for 
eigenvalues outside the extrema of SpcCggg(A). However the situation changes if 
we want to approximate an eigenvalue in a gap of SpeCgsg(A). There is no easy 
minimax principle, and spectral pollution may happen at any point of the gap. 

The difficulties involved in the computation of these eigenvalues are well 
known for particular operators, see e.g. [BoBr] or [RaSa^Va]. Moreover, in a 
generic situation we have 

Lemma 2.2. IfX ^ SpeCggg(A) is such that a < \ < (3 where a, (3 E SpeCess(A), 
there exists a sequence of subspaces Cn satisfying (2.2) for all p E N and all 
u E Dom(A), such that A E Spec{AcJ for all n E N. 

This lemma directly follows from [LeSh, Theorem 2.1]. 

2.2 The abstract quadratic projection method 

Let, as before, £ be a finite-dimensional subspace of Dom{A), and let E — 
{ei, . . . , Cn} be a basis of C. This basis need not be orthogonal. 
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Consider the quadratic matrix polynomial 

Pc{z) ■.= Qc-2zAc + z^Bc, 



(2.3) 



where 

[Bc\jk = {ej, Ck) [Ac]ik = {Acj, Ck) [Qc]jk = {Acj, Ack) . (2.4) 

In numerical analysis, Ac is called the stiffness matrix, Be is a mass matrix, and 

Qc is a bending matrix. If E is an orthonormal basis, then Ac = H^A \ C, 
and Bc = ld \ C. Additionally, if E C Dom(A2), then Qc = n^^^ ^ £ gnd 

p^(z) = UciA - r ^. 

We define the spectrum of the matrix polynomial Pc, Spec(Pc), as the set 
of // e C such that 

Pc{^i)u = for some ueC\ {0}. (2.5) 

Since Be is non-singular, det{Pciz)) is a polynomial in z of degree 2dim(£). 
Moreover, if /i G Spec(P£), then also fl E Spec(P£). Therefore Spec(P/;) is a 
set of at most 2dim(£) complex points, symmetric with respect to the real axis. 

The core idea of the quadratic projection method lies in the fact that Spec(74) 
can be well estimated if one knows the points of Spec(P£) which are "close" 
to the real line, see Corollary 2.5 and Theorem 2.7 below. In [Sh], [LeSh] and 
[Bol], Spec(Pc) is called the second order spectrum of A relative to C. This set 
was first studied in connection with the spectrum of A in [Da], where the name 
originated. 

Remark 2.3. Intuitively, the quadratic projection method arises from the follow- 
ing simple observation. Let C £ R lie in a gap of the essential spectrum. By 
virtue of the spectral theorem, the discrete eigenvalues of {( — A) inside the 
corresponding shifted gap of {( — A) containing the origin, are also the discrete 
eigenvalues of {( — A)'^ lying below the bottom of the essential spectrum of 
(C — Ay. This suggests that the truncations of the latter operator must pro- 
vide information about the localisation of a portion of SpeCfjigj,(A) near The 
quadratic projection method is a rigorous realisation of a similar idea. 

The main reason for preferring (2.5) over (2.1) for estimating the spectrum of 
A lies in the following observation. Let D{a,b) be the open disk in the complex 
plane with an interval [a,b] as a diameter: 

6 — a 1 



|w e C : 


a + b 


w 




2 
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Theorem 2.4 ([LeSh, Lemma 5.2]). Suppose that {a,b) n Spec(A) = 0. If 
z e D{a, h), then the matrix Pc{z) is non-singular. 

Proof. Our proof is slightly different from that of [LeSh]. Let z ^ D{a, h). Let 

■.= {{X- zf : A e (-00, a] U [6, oo)}. 

We first show that ^ ConvE^ (here ConvQ denotes the convex hull of the 
set Q. C C). Indeed, let 9 be the angle at z of the triangle T whose vertexes 
are a. h. z. Elementary geometric arguments show that 6 > n/2. Then the 
transformation m : A i— > (A — z)'^, maps the angular region 

B — {{w — z) : pw E T for some p > 0} 

into another angular sector centred at the origin with angle 29 > %. Since 

(-00, a) U (6, oo) C C \ S and 

m : (— oo, a] U [b, oo) i — > E^, 

there exists — tt < 9o < n and c > 0, such that Re(e'^°'u;) > c for all w e E^. 
This ensures that ^ ConvE^ as required. 

Since A = A*, {A — z)^ with domain T)om{^A?) is a normal operator, [Ka]. 
As we have for the numerical range 

Num {A - zf C Conv [Spec(^ - zf] C Conv E^ , 

and Dom(/l^) is a core for A, we have 

Re [e'^°{{A - z)u, {A - z)u)] > c 

for all u e Dom(y4) with = 1. In particular this holds true for a\\ u E C 
with \\u\\ = 1, so that Pciz) cannot be a singular matrix. □ 

As a consequence of Theorem 2.4, the points of Spec(Pc) which are close to 
the real line, are necessarily close to Spec(yl). In other words, the method never 
pollutes. We also have two immediate corollaries. 

Corollary 2.5. If e Spec(Pc), then 

inf{|Re/i-A| : A e Spec(A)} < | Im//|. (2.6) 
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If A G Spec(A) is isolated from other point of the spectrum, (2.5) provides a 
two-sided estimate of A, with an error explicitly determined without the need for 
computing eigenfunctions. In case this error is small, we can actually improve it 
by a square: 

Corollary 2.6. Let A G Spec(>l). Assume that A is isolated from other points 
of the spectrum and let 



S :— min{|A — u\ : v ^ Spec(74)} 

= dist(A,Spec(A)\{A}). 



(2.7) 



If <5/2 for /X e Spec(Pc), then 



Re/x- A| < , ' . (2.8) 



Proof. Theorem 2.4 yields 

5 



>\. (2.9) 



Using the assumption |/x — A| < 5/2, we can re-write (2.9) as 



Thus 



5 

Re/x-A| < - -y--(Im/x)2. 



Re/i — A < , < ; 

2 + Y 



□ 



Corollary 2.6 supersedes Corollary 2.5 once we have found points of Spec(Pc) 
sufficiently close to an isolated point of the spectrum of A. Note that A e 
Spec(yl) does not have to be a discrete eigenvalue. 

The above "non-pollution" results are useful as long as there are points of 
Spec(Pc) near to the real line. It is not immediately clear, however, whether or 
not the eigenvalues of A are approximated by some points in Spec(P£) when the 
dimension of C goes to infinity. The results of [Bol] and [Bo2] show that this is 
indeed the case, under a condition analogous to (2.2). 
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Theorem 2.7 ([Bo2, Theorem 2.2]). Let A G Spcc^i.c!^), and let Sx := {u : 

Au = Xu} be the corresponding eigenspace. Let L„ C Dom(A^) he suhspaces 
with corresponding orthogonal projections TLc„, such that (2.2) holds for p = 
0, 1, 2 and all u e 8\. Then there exist eigenvalues A„ e Spec(P£„) such that 
A„ — > A as n — > oo. 

3 The quadratic projection method for perturbed 
periodic Schrodinger operators 

Let H be the differential expression defined by (11) acting on the dense domain 
Let 

p = 2 if < 3, 
p>2 if = 4, 
p> N/2 if AT > 5. 

Below and elsewhere we assume that the potential V : M.^ — > M is uniformly 
locally LP in the sense that 

/ \V{x)fd^x<M (3.1) 
Jc 

for any unit hyper-cube C, where the constant M is independent of C. 

The condition (3.1) ensures that the operator of multiplication by V is (— A)- 
bounded with relative bound equal to 0, so that H '\s a self-adjoint operator and 
C^(M^) is a core for H (cf. [ReSi, Theorem XIII. 96]). Furthermore, H is 
bounded below. 

3.1 Approximating subspaces in the quadratic projection 
method for the Schrodinger operator 

We have already established, in the abstract setting of Theorem 2.4 that, for any 
choice of a subspace £ C VF^'^(R^), the eigenvalues of the matrix polynomial 
Pciz) lying close to the real axis will be close to the spectrum of H (and those 
"far away" from the real axis don't matter). In other words, the quadratic pro- 
jection method does not pollute. In order, however, to achieve a small error and 
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approximate as many eigenvalues as possible, the choice of £ (or of a sequence 
of such spaces) is absolutely crucial, see Theorem 2.7. Two main difficulties here 
are the infinite geometry and the extra smoothness requirements needed for Qc 
to make sense, see (2.4). 

Let := [— Let M.^ '■— Wq''^{Qs) be a nested family of Sobolev 
spaces. Let n G N, be a sequence of n-dimensional subspaces of Ms- Let 

{(l>s,n,k}k=i be a basis for Ms,n- Set, for j, A; = 1, . . . , n, 

[-^s,n]j,fe ■ / 4's,n,j4's,n,k i 

[As,n]j,k ■= / ^4>s,n,j ■ ^(l)s,n,k + V (ps,n,j<Ps,n,k , (3.2) 

,n,j4's,n,k ■ 

and consider a quadratic {n x n)-matrix polynomial 

-fs,n(^) Qs,n "^^-^SjU ~l~ -^s,n • ('^■^) 

Now, let Sn be a monotone increasing unbounded sequence of positive real 
numbers, let £„ = Ms„.n, and let Pn{z) = Ps„,n{z). Then Theorem 2.7 still 
holds as long as one can verify (2.2) for p = 0,1, 2. 

If the potential V is sufficiently smooth, a natural choice of the basis func- 
tions 4>s,n,k afs piecewise splines on fig satisfying (plan^ — d4)/dn\Qa^ ~ 0. 
However, even for this simple choice, verifying (2.2) is still highly technical, and 
we omit the details. 

Even fixing both parameters s and n and not imposing any condition on M.s,n 
except M.s,n C Wq''^{VLs), still usually provides some useful information about 
the spectrum, with a posteriori two-sided estimates: if A„ e Spec(Ps,n) and 
dist(Re A„, SpeCggg(if)) > | ImA„|, then there exists A e SpeC(jigp(if) which lies 
in the the same spectral gap as RcA„. See Corollary 2.6 for a sharper estimate. 

On the other hand, to achieve approximation it is crucial that both parameters 
n and s go to infinity in our choice of approximate spaces If we fix an 
arbitrarily large s and let n — > oo, then, though we still do not have pollution 
(unlike a standard projection method), neither we have approximation. 
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3.2 The quadratic matrix polynomial problem 



The quadratic projection method prescribes finding the spectrum Spec(P) of the 
a quadratic matrix polynomial of the form 

P(z) = Q + 2zA + z^B, 

cf. Section 2.2. In applications, the matrix coefficients Q, A and B are expected 
to be sparse and real. They are always hermitean, so P{z) is a self-adjoint matrix 
polynomial in the sense of [Go]. 

The standard way of finding Spec(P) is to construct a suitable companion 
linear pencil eigenvalue problem, 

Lv = iiKv for some v e C, (3.4) 

such that iJ, e Spec(P) if and only if (3.4) holds true. The coefficients, L, K, 
of the companion form, (L — zK), are twice the size of the coefficients of P{z). 
They are not unique. Two possible companion forms are given by: 

where N \s a non-singular matrix. 

Different companion forms lead to different stability properties of the linear 
pencil problem to be solved once the matrices have been assembled. It is desirable 
finding a companion form that does not worsen the condition numbers of the 
original matrix polynomial spectral problem. For a thorough account on this issue 
see [HiMaTi] and references therein. 

3.3 Examples 

One-dimensional example — Gaussian perturbation of the Mathieu po- 
tential 

Let 

= 1, Vp{x) = cos(x), Vd{x) = -e-^', 
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and H as in (1.1) with potential V = Vp + Vd- We now illustrate how to 
implement the theoretical discussion carried out in the previous sections to the 
study of Spec{H). 

The essential spectrum of H is determined by Vp. It comprises an infinite 
number of non-intersecting bands [q;„,/3„] of absolutely continuous spectrum 
whose endpoints are determined by the Mathieu characteristic values [In, §7.4]. 
The approximate endpoints of the first five bands are given in Table 1. 



n 




Pn 


1 


-0.378490 


-0.347670 


2 


0.594800 


0.918058 


3 


1.29317 


2.28516 


4 


2.34258 


4.03192 


5 


4.03530 


6.27082 



Table 1: Endpoints of the first five bands of the essential spectrum for the 
Gaussian perturbation of the Mathieu potential. Cf. [AbSt]. 

Addition of the negative Gaussian potential yields a non-empty discrete spec- 
trum. By implementing the quadratic projection method (3.3) into a finite ele- 
ment scheme, we detect three eigenvalues of H with high accuracy: 

Ai ^ -0.40961, A2 ^ 0.37763, A3 ^ 1.18216. 

The eigenvalue Ai is below the bottom of the essential spectrum, whereas A2 
and A3 lie in the first and the second gap, respectively. 

Figure 1 illustrates the main ideas discussed in the previous sections. The 
spectrum of P{z) is shown as blue dots, while the eigenvalues of the standard 
Galerkin eigenvalue problem (2.1) are shown as red crosses. The picture shows 
a narrow strip of the complex plane with the bottom edge being the interval 
[—0.5, 2]. Note that there are eigenvalues of P{z) close to each of the eigenvalues 
Ai, A2 and A3. According to Corollary 2.5, these eigenvalues are not spurious: the 
real part of a complex number z G Spec(P) is always an approximation of points 
in Spcc(i/), with a two-sided error estimate depending on Imz. There are also 
eigenvalues of the linear problem (2.1) near to Spec^ig^{H). These eigenvalues 
also provide one-sided approximation (from above) for the Xj. However one 
should be careful when using the Galerkin method, as spectral pollution may 
happen. For this particular set of parameters there are two spurious eigenvalues: 
one near —0.3 and the other near 1.1. 
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Figure 1: The quadratic projection method vs the Galerkin projection method. 
Here s = 49. Insert: zoom near A ^ —0.35. 

Two-dimensional examples 

We now consider a family of case studies with N = 2. For / G W^^'^(R^), let 

Hof{x,y) = -Af{x,y) + (cos(x) + cos{y))f{x,y) 
HJ{x, y) = Hofix, y) - ce-(^'+^')/(x, y) 
H^f{x,y) = Hof{x,y) - cxe-^^'^y'^ f{x,y), 

where c > 0. A straightforward argument involving separation of variables shows 
that 

Spec(iJo) = [J {A + /i : /i G Spec(if)}, 

AGSpec(_ff) 

where H = —d'^ + cos{x) is the one dimensional Mathieu Hamiltonian. Further- 
more, as both Hi and H2 are relatively compact perturbations of Hq, 

SpeCess(^i) = Spec^33(if2) = Spec{Ho). 
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An approximation of the endpoints of the bands comprising the essential spec- 
trum is given in Table 2. Unlike the one-dimensional model, we now have a 
finite number of gaps. Note that the perturbation associated to Hi is radially 
symmetric and sign definite, while the one associated to H2 is sign indefinite and 
not radially symmetric. 



n 


On 




1 


-0.756978 


-0.695338 


2 


0.216310 


0.570389 


3 


0.914677 


00 



Table 2: Endpoints of SpeCess(i7o)- 



With the quadratic projection method we have been able to detect some 
discrete eigenvalues of Hi and H2 for different values of the coupling constant c. 
These results are presented in Table 3. As we increase the value of c, eigenvalues 
of Hi are moving from right to left. From the numerical results, the same seems 
to be true for eigenvalues of H2. Note that if an eigenvalue is close to an end- 
point of a band of the essential spectrum, the estimate (2.6) does not allow us 
to distinguish between this eigenvalue and the end-point of the band — thus the 
gaps in Table 3. 





Eigenvalues 


of Hi 








Eigenvalues 


of i?2 




c 


Ai 


A2 


C 


Ai 


A2 


5.0 




-0.09697 ± 3.39 


10" 


-4 


10 


0.1377 ±6.61 • 10"^ 




0.7559 ± 7.07 


10"^ 


5.2 




-0.17133 ±2.03 


10" 


-4 


11 


0.0865 ± 2.79 • lO'^ 




0.681 ± 1.05 


10-1 


5.4 




-0.25255 ± 1.45 


10" 


-4 


12 


0.0115 ±1.28 -10-3 








5.6 




-0.33905 ± 1.44 


10" 


-4 


13 


-0.09190 ±8.19 • 10"^ 








5.8 




-0.42902 ± 1.71 


10" 


-4 


14 


-0.22250 ± 7.77 • 10~^ 








6.0 


-0.76946 ± 2.85 • lO'^ 


-0.51978 ±2.52 


10" 


-4 


15 


-0.3730 ±1.17- 10-3 








6.2 


-0.78612 ± 8.43 • 10"^ 


-0.60527 ± 5.21 


10- 


-4 


16 


-0.5279 ± 1.97 • 10-3 









Table 3: Approximated eigenvalues of Hi and H2. 



Note that an eigenvalue Ai of the Hamiltonian Hi is below the bottom of 
the essential spectrum for c ^ 6. The Galerkin method could actually be imple- 
mented to approximate this eigenvalue. The quadratic projection, however, works 
whether an eigenvalue is in a gap or not, and also provides a good approximation 
in this case. 
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In Figure 2 we show the portion of of Spec(P) lying in the box [—1, 1/2] x 
[—3/2, 3/2] for H = H2, c = 14 and s = 60. Corresponding pictures for Hi and 
other choices of c and s are qualitatively similar. This graph clearly indicates 
approximation to an eigenvalue Ai ~ —0.2225 (see the right hand picture). A 
large portion of Spec(P) forms an annular cloud around the spectral gap (Pi, 02) 
and is sufficiently away from M to indicate that there are no other eigenvalues in 
this gap. Note also that some eigenvalue of P{z) are close to SpeCggg(if2)- 



0.5 



-0.5 



, 7/ • , ...... _ 



-0.5 



O.J 



0.09- 
0.08- 
0.07- 
0.06- 
0.05- 
0.04- 
0.03- 
0.02- 
0.01 - 
0- 

-0.01 - 



-0.5 



0.5 



Figure 2: The quadratic projection method for our two dimensional models. Left: 
Typical output in the computation of Spec(P) for operator H2 when c = 14 (here 
s = 60). Right: zoom in the left picture on a narrow strip near the real line. 



4 Final remarks 

Other procedures exist for computing the eigenvalues of perturbed periodic partial 
differential operators such as H, see [Do]. These include a method based on 
finding the eigenvalues of the matrix pencil problem 

As,nU = XBs^nU foY some M G £ \ {0}, (4.1) 

where the matrix coefficients are defined by (3.2) (that is applying the projection 
method) for several values of s and n, and observing the dynamics of the eigen- 
values of (4.1) as s increases. Some of the eigenvalues of (4.1) will be spurious 
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and some will be close to the true spectrum of H. The spurious eigenvalues will 
typically be unstable as functions of the parameter ,s. The approximate eigen- 
values close to the true spectrum of H will be, on the other hand, very stable. 
Thus, by increasing s, and tracking the evolution of the eigenvalues of (4.1), one 
would be able to obtain some information about Spec(if). 

This method, however, is quite inaccurate and it becomes useless when N > 
2, and we are interested in finding large eigenvalues. Furthermore, it very much 
depends upon the choice of approximating subspaces We are not aware of 
any rigorous treatment of the effectiveness of this approach. 

As the chosen subspaces Wq''^{Qs) are naturally nested for increasing values 
of s > and they are all embedded in l^^'^(R^), every point in Spec(if) is 

approximated (always from above) by the spectrum of (H \ Wq''^{Qs))- Note 
that compactly supported functions form a core for the operator, and satisfy any 
boundary condition if the boundary is far enough away. Spectral pollution in 
the projection method is a consequence of high eigenvalues of {H \ Wq''^{Qs)) 
accumulating at the bottom of the essential spectrum of H, and this effect is 
unavoidable. 

We suggest using instead (or in addition to standard techniques), the quadratic 
projection method, which never pollutes. 

Acknowledgements 

We are grateful to Marco Marietta for useful discussions and helpful advice. 

References 

[Bi] M. BiRMAN, The discrete spectrum in gaps of the perturbed periodic 

Schrodinger operator. I. Regular perturbations. Boundary value prob- 
lems, Schrddinger operators, deformation quantization, Math. Top., 
Akademie Verlag, Berlin, 8 (1995) 334-352. 

[AbSt] M. Abramowitz, I. Stegun, Handbook of mathematical func- 
tions with formulas, graphs and mathematical tables, National Bureau 
of Standards, 1964. 

[BoBr] D. BOFFi, F. Brezzi, L. Gastaldi, On the problem of spurious 
eigenvalues in the approximation of linear elliptic problems in mixed 
form. Math. Comp. 69 (1999) 121-140. 



15 



[Bol] L. BouLTON, Limiting set of second order spectrum. Math. Comp. 
75 (2006) 1367-1382. 

[Bo2] L. BouLTON, Non-variational approximation of discrete eigenvalues 
of self-adjoint operators, to appear in IMA J. Numer. Anal. (2006). 

[Da] E.B. Davies, Spectral enclosures and complex resonances for gen- 

eral self-adjoint operators. LMS J. Comput. Math. 1 (1998) 42-74. 

[DeHe] P. A. Deift, R. Hempel, On the existence of eigenvalues of the 
Schrodinger operator H — XW in a gap of cr(if). Commun. Math. 
Phys. 103 (1986) 461-490. 

[DoEsSe] J. DoLBEAULT, M.J. EsTEBAN, E. Sere, On the eigenvalues of 
operators with gaps. Application to Dirac operators. J. Funct. Anal. 
174 (2000) 208-226. 

[Do] J. P. DOWLING, Photonic &. Sonic Band-Gap Bibliography, 

http : //phys . Isu . edu/~ jdowling/pbgbib . html 

[Go] I. GoHBERG, p. Lancaster, L. Rodman, Matrix Polynomials, 
Academic Press, New York, 1982. 

[HiMaTi] N. Higham, S. Mackey, F. Tisseur The Conditioning of Lin- 
earizations of Matrix Polynomials. SI AM J. Matrix Anal. Appl. (2006) 
to appear. 

[In] E. Inge, Ordinary Differential Equations, Dover, New York, 1956. 

[Ka] T. Kato, Perturbation theory for linear operators, 2nd edition, 

Springer-Verlag, Berlin, 1980. 

[Ku] P. KuGHMENT, Floquet Theory For Partial Differential Equations, 

Birkhauser Verlag, Basel, 1993. 

[LeSh] M. Levitin, E. Shargorodsky, Spectral pollution and second 
order relative spectra for self-adjoint operators. IMA J. Numer. Anal. 
24 (2004) 393-416. 

[RaSa^Va] J. Rappaz, J. Sanghez Hubert, E. Sanchez Palengia, D. 

Vassiliev, On spectral pollution in the finite element approximation 
of thin elastic 'membrane' shell. Numer. Math. 75 (1997) 473-500. 



16 



[ReSi] M. Reed, B. Simon, Methods of Modern Mathematical Physics, 
Volume 4: Analysis of Operators, Academic Press, New York, 1978. 



[Sh] E. Shargorodsky, Geometry of higher order relative spectra and 

projection methods. J. Oper. Theo. 44 (2000) 43-62. 

[StFi] G. Strang, G. Fix, An Analysis of the Finite Element Method, 
Prentice-Hall, New Jersey, 1973. 

[Su] T. SuSLiNA, The discrete spectrum of a two-dimensional second- 

order periodic elliptic operator perturbed by a decaying potential. II. 
Inner gaps. Algebra i Analiz 15 (2003) 128-189; translation in St. 
Petersburg Math. J. 15 (2004) 249-287. 



17 



